home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 5 / Apprentice-Release5.iso / Source Code / Libraries / Advanced I⁄O v2.3 / Advanced i⁄o / arithm_modadh.cc < prev    next >
Text File  |  1995-07-08  |  10KB  |  328 lines
  1. // This may look like C code, but it is really -*- C++ -*-
  2. /*
  3.  ************************************************************************
  4.  *
  5.  *                    Arithmetic Coding
  6.  *           Adaptive Histogram-based model for the data source
  7.  *
  8.  *
  9.  * The present file provides the functionality for both building the
  10.  * histogram for a sequence of integers items and supplying the
  11.  * arithmetic codec with probabilities of the items estimated from
  12.  * the histogram.
  13.  *
  14.  * The program uses the histogram to find out which symbols shall really
  15.  * come out and how frequent. It allocates frequency tables only for
  16.  * the symbols that really rather than potentially need to be coded.
  17.  * In all other respects (assuming Lorentzian-like initial probability
  18.  * distribution and tailoring as more symbols become known) the
  19.  * model behaves way the same as the simple AdaptiveModel does.
  20.  *
  21.  * $Id: arithm_modadh.cc,v 2.0 1995/02/07 19:37:03 oleg Exp oleg $
  22.  *
  23.  ************************************************************************
  24.  */
  25.  
  26. #ifdef __GNUC__
  27. #pragma implementation "arithm_modadh.h"
  28. #endif
  29.  
  30. #include "arithm_modadh.h"
  31. #include "std.h"
  32.  
  33.  
  34. /*
  35.  *========================================================================
  36.  *        Histogram-based Adaptive model for the data source
  37.  */
  38.  
  39.  
  40. /*
  41.  *------------------------------------------------------------------------
  42.  *            Constructors and Destructors
  43.  */
  44.  
  45.                 // Constructor
  46. AdaptiveHistModel::AdaptiveHistModel(const Histogram& histogram)
  47.      : no_distinct_symbols(histogram.no_distinct_symbols()),
  48.        update_inc(0)
  49. {
  50.                  // Skip leading and tailing dummy symbols
  51.                  // (never occured in the histogram)
  52.   for(symbol_lwb=histogram.symbol_lwb;
  53.       symbol_lwb < histogram.symbol_upb && histogram.get(symbol_lwb) == 0;
  54.       symbol_lwb++)
  55.      ;
  56.      
  57.   for(symbol_upb=histogram.symbol_upb;
  58.       symbol_upb > symbol_lwb && histogram.get(symbol_upb) == 0;
  59.       symbol_upb--)
  60.      ;
  61.   assure( symbol_upb > symbol_lwb,
  62.         "there should be at least 2 symbols to code" );
  63.   no_potential_symbols = symbol_upb - symbol_lwb + 1;
  64.   assert( no_potential_symbols >= no_distinct_symbols );
  65.   allocate_model(no_distinct_symbols);
  66.   symbol_to_index = new Index[no_potential_symbols];
  67.   index_to_symbol = new Symbol[no_distinct_symbols];
  68.  
  69.   initialize_model(histogram);
  70. }
  71.  
  72.                 // Constructor for a dummy model. The real
  73.                 // model contents is supposed to be read
  74.                 // in later
  75. AdaptiveHistModel::AdaptiveHistModel(void)
  76.      : symbol_lwb(0), symbol_upb(0),
  77.        no_potential_symbols(0),
  78.        no_distinct_symbols(0),
  79.        update_inc(0)
  80. {
  81.   symbol_to_index = 0;
  82.   index_to_symbol = 0;
  83. }
  84.  
  85.                 // Destructor
  86. AdaptiveHistModel::~AdaptiveHistModel(void)
  87. {
  88.   assert( symbol_to_index != 0 && index_to_symbol != 0 );
  89.   delete [] symbol_to_index;
  90.   delete [] index_to_symbol;
  91. }
  92.  
  93. /*
  94.  *------------------------------------------------------------------------
  95.  *           Sort the histogram and set up the model
  96.  */
  97.  
  98. static struct Hist_elem 
  99.     { unsigned int count; Symbol symbol; }
  100. * Hist_sorted;
  101.  
  102.                 // Comparison routine to sort the Histogram
  103.                 // in the descending order
  104. static int hist_compar(const void * eli, const void * elj)
  105. {
  106.   if( ((const Hist_elem*)eli)->count < ((const Hist_elem*)elj)->count )
  107.     return 1;
  108.   else if( ((const Hist_elem*)eli)->count == ((const Hist_elem*)elj)->count )
  109.     return 0;
  110.   else
  111.     return -1;
  112. }
  113.  
  114. void AdaptiveHistModel::initialize_model(const Histogram& histogram)
  115. {
  116.   Hist_sorted = new Hist_elem[no_distinct_symbols];
  117.  
  118.   register Symbol symbol;
  119.   register int i;
  120.   for(symbol=histogram.symbol_lwb,i=0; symbol<=histogram.symbol_upb; symbol++)
  121.   {
  122.     unsigned int count = histogram.get(symbol);
  123.     if( count != 0 )
  124.       Hist_sorted[i].count = count, Hist_sorted[i].symbol = symbol, i++;
  125.   }
  126.   assert( i == no_distinct_symbols );
  127.  
  128.   qsort(Hist_sorted,no_distinct_symbols,sizeof(Hist_sorted[0]),hist_compar);
  129.  
  130.   memset(symbol_to_index,0,sizeof(symbol_to_index[0])*no_potential_symbols);
  131.   for(i=0; i<no_distinct_symbols; i++)
  132.   {
  133.     Index  index  = i+1;
  134.     Symbol symbol = Hist_sorted[i].symbol;
  135.     assert( symbol <= symbol_upb && symbol >= symbol_lwb );
  136.     symbol_to_index[symbol-symbol_lwb]  = index;
  137.     index_to_symbol[index-1] = symbol;
  138.   }
  139.   assert( i<EOF_index );
  140.  
  141.   initial_distribution();
  142. }
  143.  
  144.                 // Assign initial frequency counts
  145.                 // according to the Lorentzian ditribution
  146. void AdaptiveHistModel::initial_distribution(void)
  147. {
  148.   register int i;
  149.   cum_frequencies[no_indices] = 0;
  150.                     // update_inc > 1 means the real
  151.                     // occurence of a symbol means more
  152.   update_inc = no_indices/10 + 1;     // than just a potential occurence.
  153.                     // Since all frequencies must be >0,
  154.                     // frequency count one can be regarded
  155.                     // as being "infinitely small"
  156.  
  157.                     // The folowing distribution is
  158.                     // certaily ad hoc, but seems to work
  159.   for(i=no_indices; i>0; i--)
  160.   {
  161.     frequencies[i] = ( i <= 2 ? 45*update_inc-30 : 
  162.                    i <= 12 ? 3*update_inc : 1 );
  163.     cum_frequencies[i-1] = cum_frequencies[i] + frequencies[i];
  164.   }
  165.   frequencies[0] = 0;            // Must be zero (an impossible value)
  166.  
  167.  
  168. }
  169.  
  170. /*
  171.  *------------------------------------------------------------------------
  172.  *              Searching for index/symbol
  173.  */
  174.  
  175.                 // Return the index of a given symbol
  176. Index AdaptiveHistModel::get_index(const Symbol symbol) const
  177. {
  178.   if( symbol < symbol_lwb || symbol > symbol_upb )
  179.     _error("Symbol %d is out of interval [%d,%d]",
  180.        symbol, symbol_lwb, symbol_upb);
  181.   Index index = symbol_to_index[symbol-symbol_lwb];
  182.   if( index == 0 )
  183.     _error("\nSymbol %d has not been counted when Histogram was constructed",
  184.        symbol);
  185.   return index;
  186. }
  187.  
  188.                 // Return the symbol for a given index
  189. Symbol AdaptiveHistModel::get_symbol(const Index index) const
  190. {
  191.   if( index < 1 || index > no_indices || index == EOF_index )
  192.     _error("Index %d is out of boundaries [1,%d)",index,no_indices);
  193.   return index_to_symbol[index-1];
  194. }
  195.  
  196. /*
  197.  *------------------------------------------------------------------------
  198.  *           Input/Output for the frequency tables
  199.  * As a matter of fact, only  index_to_symbol[] table along with its size 
  200.  * needs to be saved. All other tables can easily be reconstructed from that.
  201.  *
  202.  * Thus, the format the tables are written is follows
  203.  *    short    no_distinct_symbols
  204.  *    short    index_to_symbol[0..no_distinct_symbols-1]
  205.  * The tables entries are coded using a variable bit size algorithm,
  206.  * see BitIn/BitOut for more details.
  207.  */
  208.  
  209. void AdaptiveHistModel::open(BitOut& file)
  210. {
  211.   file.write_short(no_distinct_symbols);
  212.   for(register int i=0; i<no_distinct_symbols; i++)
  213.     file.put_short(index_to_symbol[i]);
  214. }
  215.  
  216.                     // Read in the tables and set up
  217.                     // the model
  218. void AdaptiveHistModel::open(BitIn& file)
  219. {
  220.   no_distinct_symbols = file.read_short("Reading AdaptiveHistModel tables");
  221.   allocate_model(no_distinct_symbols);
  222.   index_to_symbol = new Symbol[no_distinct_symbols];
  223.  
  224.   for(register int i=0; i<no_distinct_symbols; i++)
  225.   {
  226.     const Symbol symbol = file.get_short();
  227.     if( i== 0 )
  228.       symbol_lwb = symbol_upb = symbol;
  229.     else
  230.       symbol_lwb = symbol_lwb > symbol ? symbol : symbol_lwb,
  231.       symbol_upb = symbol_upb < symbol ? symbol : symbol_upb;
  232.     index_to_symbol[i] = symbol;
  233.   }
  234.   
  235.   no_potential_symbols = symbol_upb - symbol_lwb + 1;
  236.   assert( no_potential_symbols > 1 );
  237.   assert( no_distinct_symbols <= no_potential_symbols );
  238.  
  239.   symbol_to_index = new Index[no_potential_symbols];
  240.  
  241.   memset(symbol_to_index,0,sizeof(symbol_to_index[0])*no_potential_symbols);
  242.   for(register Index index=1; index<=no_distinct_symbols; index++)
  243.   {
  244.     Symbol symbol = index_to_symbol[index-1];
  245.     symbol_to_index[symbol-symbol_lwb] = index;
  246.   }
  247.   
  248.   initial_distribution();
  249. }
  250.  
  251. /*
  252.  *------------------------------------------------------------------------
  253.  *            Update the adaptive model
  254.  *        to account for occurence of a particular index
  255.  */
  256.  
  257.                 // Scale the accumulated statistics down
  258.                 // in anticipation of a change, or simply to
  259.                 // keep counters from overflow
  260. void AdaptiveHistModel::scale_down_past(void)
  261. {
  262.   register int cum = 0;            // If so, halve all the counts
  263.   for(register int i=no_indices; i>=0; i--)
  264.   {                    // keeping them nonzero, and
  265.     cum_frequencies[i] = cum;        // recompute the cumulative counts
  266.     cum += (frequencies[i] = (frequencies[i]+1) >> 1);
  267.   }
  268. }
  269.  
  270. void AdaptiveHistModel::update_model(const Index index)
  271. {
  272.                     // See if freq counts near their max
  273.   if( total_cum_freq() >= (Top_freq>>1) ) // /2 makes rescaling work more often
  274.     scale_down_past();                    
  275.                       // Tentatively increment the freq
  276.                       // count for index, just an estimate
  277.   const Lword tent_freq = frequencies[index] + 
  278.                   (update_inc > 4 ? update_inc : 1);
  279.  
  280.                     // Try to pull 'index' to the beginning
  281.                     // of the freq. table, skipping symbols
  282.                     // with the lower freq.
  283.   register Lword * fp = &frequencies[index];
  284.   while( *--fp != 0 && *fp <= tent_freq )    // Note, freq[0] is always 0
  285.     ;
  286.                       // Find a new position for the index
  287.   const Index new_index = ++fp - frequencies;
  288.   if( new_index < index )        // If the index is to move, update
  289.   {                    // the translation tables
  290.     Symbol symbol_oind = index_to_symbol[index-1];
  291.     Symbol symbol_nind = index_to_symbol[new_index-1];
  292.     index_to_symbol[index-1] = symbol_nind;
  293.     index_to_symbol[new_index-1] = symbol_oind;
  294.     symbol_to_index[symbol_oind-symbol_lwb] = new_index;
  295.     symbol_to_index[symbol_nind-symbol_lwb] = index;
  296.  
  297.                         // Note, we've just swapped index
  298.                         // and new_index. If their frequencies
  299.                         // differ, we need to update the
  300.                         // frequecies and cumulative freqs
  301.     const Lword old_freq = frequencies[index];
  302.     const Lword new_freq = frequencies[new_index];
  303.     if( old_freq != new_freq )
  304.     {
  305.       frequencies[index] = new_freq;
  306.       frequencies[new_index] = old_freq;
  307.       const long int delta = new_freq - old_freq;
  308.       for(register Lword * cfp = &cum_frequencies[new_index];
  309.           cfp < &cum_frequencies[index];)
  310.          *cfp++ += delta;
  311.     }
  312.   }
  313.  
  314.   frequencies[new_index] += update_inc;    // Increment the freq count for index
  315.   register Lword * cfp;            // And update the cumulative counts
  316.   for(cfp=&cum_frequencies[new_index]; cfp>&cum_frequencies[0]; )
  317.     *--cfp += update_inc;
  318.  
  319. #if 0
  320. {
  321.   message("index %d, new_index %d, update_inc %d\n",index,new_index,update_inc);
  322.   for(register Index i=0; i<no_indices; i++)
  323.      message("%d ",frequencies[i]);
  324.   message("\n");
  325. }
  326. #endif
  327. }
  328.